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ABSTRACT 

We  discuss  two  methods  for  reducing  the  variance  of 
estimators  of  parameters  of  limiting  distributions  of  stable 
stochastic  processes  in  simulations.   The  methods  are  discussed 
in  the  context  of  the  simple  GI/G/1  queue.   Of  the  two  methods 
one,  which  we  call  an  internal  control  variable,  gives  a  vari- 
ance reduction  which  is  roughly  uniform  over  values  of  the 
parameters  of  the  process  and,  in  particular,  works  well  for 
values  of   p.   the  traffic  intensity,  close  to  1. 
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VARIANCE  REDUCTION  FOR  REGENERATIVE  SIMULATIONS,  I: 
INTERNAL  CONTROL  AND  STRATIFIED  SAMPLING  FOR  QUEUES 

by 

Donald  L.  Iglehart      and      Peter  A.  W.  Lewis 

1.   Introduction  and  Summary 

Simulators  are  frequently  faced  with  the  task  of  estimating 
a  parameter  associated  with  the  limiting  distribution  of  a  stochas- 
tic process  which  is  being  simulated.   A  methodology  based  on 
regenerative  processes  for  obtaining  point  estimates  and  confidence 
intervals  for  such  parameters  has  recently  been  developed  in  Crane 
and  Iglehart  (1974a, b) ,  (1975a, b)  and  Iglehart  (1975) ,  (1976a, b) . 
In  this  paper  we  shall  indicate  two  techniques,  internal  control 
variables  and  internal  stratified  sampling,  which  might  be  used  in 
conjunction  with  the  regenerative  method  for  obtaining  additional 
variance  reduction  for  the  estimates.   To  illustrate  these  tech- 
niques we  shall  restrict  our  attention  in  this  paper  to  the  GI/G/1 
queue.   Other  applications  of  these  ideas  will  be  dealt  with  in 
future  publications,  as  well  as  uses  of  the  regenerative  property 
which  may  possibly  be  more  suited  for  obtaining  variance  reduction 
for  estimates  in  stable  stochastic  processes. 

In  the   GI/G/1   queue  we  assume  the  zeroth  customer  arrives 

at  time   tQ  = 0 ,   finds  a  free  server,  and  experiences  a  service 

time   v~  .   The   n —  customer  arrives  at  time   t   and  experiences 
0  n 

a  service  time   v  .   Let  the  interarrival  times   t  -  t    =  u  , 

n  n    n-1    n 

n>l.   Assume  the  two  sequences   {v  :  n>0}   and   {u  :  n>l}   each 

n  n 


consist  of  independent,  identically  distributed  (i.i.d.)  random 

variables  (r.v.'s)  and  are  themselves  independent.   Let  E{v  }= 

n 

\i~    ,   E(u  }  =  X~  ,   and   p=A/y,   where   0<A,vi<°°.   Thus   y  has 
the  interpretation  of  the  mean  service  rate  and   A   has  the  inter- 
pretation of  the  mean  interarrival  rate.   The  parameter   p   is 
called  the  traffic  intensity  and  is  the  natural  measure  of  conges- 
tion for  this  system.   We  shall  assume  that   p<  1,   a  necessary 
and  sufficient  condition  for  the  system  to  be  stable. 

While  many  characteristics  of  interest  can  be  estimated 
using  the  regenerative  method,  we  shall  restrict  our  attention  to 
the  waiting  time  of  the   n —  customer,   W    (time  from  arrival  to 
commencement  of  service) .   To  obtain  a  representation  for  the  pro- 
cess  {W  :  n>0}   let  X  =v   ,  -u   and  set  S  =0,   S  =X,  + +X  , 

n  n   n-1   n  0       n    1      n 

n>l.   The  following  well-known  recursive  relationship  exists  for 

the  W  's: 
n 

wn  =  °/   w  -lt  =  [w  +x  ^J+/   n>0. 
0        n+1     n   n+1   ' 

By  induction,  we  also  have 

W   =  max{S  -S,  :  k =  0 ,1 , . . . ,n) ,   n>0. 
n         n   k       iiit 

Using  the  strong  Markov  property  of  the  process   {S  :  n>0} 
it  can  be  shown  that  there  exists  a  sequence  of  integer-valued  r.v.'s 
{6k:  k>0>   such  that   $0  =  0,   3,  <  3k+,,   and  Wg  =  0   with  proba- 
bility  one.   In  other  words,  the  customers  numbered   3,   are  those 
lucky  fellows  who  arrive  to  find  a  free  server  and  experience  no  wait- 
ing in  the  queue.   The  fact  that  there  exists  an  infinite  number  of 
such  customers  in  the   GI/G/1   queue  is  a  direct  consequence  of  the 


assumption  that   p  <  1   and  the  strong  law  of  large  numbers.   If 
we  let   a,  =  3,  -  3k_ -,  i      k  >  1,   then   a,   represents  the  number  of 
customers  served  in  the   k —  busy  period  (b.p.)  and  they  are 
numbered   {3,  ,,  3,  , +1, .  .  .  ,  8,-1)  . 

Next  define  the  sequence   {Y  :  k>l}   by 

Yk  =     I     W 
3    pk-l 

4-Vi 

the  sum  of  the  waiting  times  in  the   k —  b.p.   Since  the  queue 

is  stable  for   p  <  1  we  have  W  =»  W   as   n-*-°°,   where   =>  denotes 

n 

weak  convergence. 

Our  goal  is  to  estimate   E{W}   by  simulation. 

It  is  known  that  E{W}  <  »   if  and  only  if  E{(X*)2}  <  °° .   We 
shall,  for  simplicity,  assume  that  E^vn^  <  °°     which  will  guarantee 
that  E{w}  <°°.   The  regenerative  simulation  method  is  based  on  the 
analytic  results  that  the  sequence   t^w01]^*  k-l}   is  independent  and 
identically  distributed  and  E{W} = E{Y  }/E{a1 } ,   the  ratio  of  two 

means.   This  suggests  using  the  ratio  of  estimates  of  E(Y1)   and 

1   n 
E(an)   to  estimate  E  (W)  .   Thus,  if  we  let  Y(n)  =-  \     Y   and 

1  n  n  k=1 

a(n)  =—  \      a  ,   where   n   now  denotes  the  number  of  cycles  observed, 

n  k=l   K 

then  a  ratio  estimate  of  E (W) ,   for  example,  is 


W(n)  =  Y(n)/a(n)  .  (1.1) 

(In  the  sequel  we  drop  the  dependence  on   n   unless  necessary.) 

Now  let   Z,  =  Y  -  E{W}a,,   k>l,   and  note  that  the  sequence 
{Z.  :  k  >  1}   is  i.i.d.  and  E{Z,}  =  0.   We  assume  the  variance  of   Z,  , 


E{z£}  =  a2  =  var{Yk>-2  cov{Yk,ak>E{w}  +  var{ak>E2 {W}    (1.2) 
is  finite  and  positive.   Then  one  can  easily  show  that  as   n  ■*■  °° 

n1/2[Y(n)/a(n)-E{W}] 


a/Eia  ) 


N(0,1),  (1.3) 


where   N(0,1)   is  a  mean  zero,  variance  one  normal  random  variable. 
This  result  yields  a  confidence  interval  for  E{w}   provided  we 
can  estimate   a/E{a.}.   A  variety  of  estimates  have  been  studied 
and  are  reported  on  in  Iglehart  (1975a).   Here  we  just  mention  two. 
The  so-called  classical  estimate  for   a/E{a,}   is  given  by 

S1  =  [Sl;L-2S12(Y/a)  +S22(Y/a)2]1/2/a, 

where   S.n   is  the  sample  variance  of  the   Y,  ' s ,   S_0   of  the   a,  '  s , 
11  r  k      22  k 


and 


S-2   the  sample  covariance  of  the   Y,  '  s   and   a,  *  s . 

The  second  estimate  of   a/E{cu}   is  the  jackknife  estimate 


which  is  defined  to  be 


L  =  [  I  (e.-e)2/(n-i)]1/2 

*     i=l    1 


-  -  -   l   n 

where   0  =  n  (Y/a)  -  (n-1)  (  \      Y./J   a.),   l<i<n,   and   9=-  \ 

is  called  the  jackknifed  estimator  of  E  (W)  .   Both   S..   and   S 


2 
are  strongly  consistent. 

The  jackknifing  technique  is  known  to  work  particularly 

well  as  a  confidence  interval  estimate  for  ratios;  for  a  large 

number  of  cycles   n   the  computational  problem  is  severe,  but  in 

that  case  the  technique  using  (1.3)  and   S   works  well.   For 

details  see  Miller  (1974)  and  Iglehart  (1975) . 


The  problem  addressed  in  this  paper  is  how  to  apply  variance 
reduction  techniques  with  the  ratio  estimator  W(n).   In  almost 
all  practical  situations,  where  in  particular  one  might  want  to 
compare  the  mean  waiting  times  of  two  different  queueing  systems 

(Iglehart  (19  76) ) ,  there  is  a  premium  on  achieving  the  minimum 
possible  variance  for  the  estimation  in  the  given  computing  time 

(number  of  cycles  allowed) .   Variance  reduction  techniques  for 
simulations  are  discussed  in  Kleijnen  (19  74)  and  Gaver  and  Thompson 

(19  73) ,  but  there  are  difficulties  in  applying  these  to  ratio 
estimates  and  regenerative  systems.   In  particular  variance  reduc- 
tion via  the  usual  control  variable  techniques  is  difficult.   A 
variation  of  this  technique,  which  we  call  internal  control  variables 
and  which  is  generally  useful  for  ratio  estimates,  is  introduced 
and  shown  to  give  considerable  variance  reduction  for  the  point  esti- 
mate W(n).   Another  technique,  internal  stratified  sampling,  is  also 
explored.   It  is  a  natural  technique  to  use  but  appears  to  be  diffi- 
cult to  use  with  ratios.   Moreover,  it  becomes  less  and  less 
effective  as   p  +  1,   while  the  internal  control  technique  holds  up 
well  for   p   close  to   1.   In  fact  the  internal  control  variables 
described  here  for  an  M/M/l   queue  give  a  variance  reduction  which 
is  fairly  uniform  for  all  values  of   y   and  ratios   p <  1 .   Better 
results  can  be  obtained  for  particular  cases  of  the  parameter   p. 

It  will  also  be  apparent  after  the  development  of  the  vari- 
ance reduction  techniques  in  the  next  sections  that  the  two  tech- 
niques for  confidence  interval  estimates  discussed  above  apply  to 
the  estimates  after  variance  reduction  has  been  applied. 


2 .   Internal  Control 

Two  main  methods  for  variance  reduction,  antithetic  variates 
and  control  variables  (Gaver  and  Thompson  (1973))  have  been  put 
forward  for  use  with  queues  in  the  non-regenerative  situation.   Of 
these,  antithetics  has  very  limited  utility  beyond  the  simple   M/M/l 
situation  in  which  it  is  patently  clear  how  to  generate  two  reali- 
zations of  samples  which  give  negatively  correlated  estimates. 
That  scheme  for  generating  negatively  correlated  estimates  does 
not  work  with  the  regenerative  method  because  the  regenerative 
blocks  in  the  original  realization  of  the  queue  and  the  antithetic 
realization  get  out  of  synchronization.   No  alternative  way  has 
been  found  to  utilize  antithetic  variates  with  the  regenerative 
technique  and  we  feel,  along  with  many  computer  scientists,  that 
the  technique  has  limited  validity  in  systems  simulation. 

The  technique  of  using  a  control  variable  and,  in  particular, 
a  regression-adjusted  control  variable  (Gaver  and  Thompson  (1973) 
p.  5  87)  is  much  more  broadly  applicable  in  systems  simulation, 
although  it  is  again  difficult  to  adapt  to  the  regenerative  situa- 
tion.  Briefly,  say  we  are  simulating  an   M/G/l   queue  to  estimate 
E (W)   with  the  non-regenerative  technique  of  averaging  the  first 
m  waiting  times  to  obtain  an  estimate  of  the  unknown   E (W) , 

m 


i  y     W.. 
m  j=l   ^ 


W   =  -      )   W.  .  (2.1) 

m   m 


The  same  random  numbers  used  to  generate  the   m  non-exponentially 
distributed  service  times  are  used  to  generate   m  exponential 
service  times  for  simulating  an   M/M/l   queue  whose  input  stream 


is  identical  to  that  of  the  simulated  M/G/l   queue.   One  would 
then  expect  that  if  the  G-distribution  is  not  too  different  from 
the  exponential  distribution,  successive  waiting  times,  say  W!, 
in  the   M/M/l   queue  realization  would  be  close  to  (positively 
correlated  with)  the  corresponding  waiting  times   W .   in  the 
M/G/l   queue.   Consequently,  the  average  of  the   W.  ' s ,  say  W  , 

will  be  positively  correlated  with  W   and  one  can  form  a  new 

r  J  m 

estimate 


Wm  =  W  +  3(W'-E(W')) .  (2.2) 

mm     mm 

Now  E(W')   is  close  to  E{W  }  =  p/y  (1-p)   for  m   large,  so 

E(W  )  ~ E (W  )   and  the  variance  of  the  new  estimate  will  be  a  minimum 

m  ~    m 


if 


s\j  f^j 


in  fact 


3  =  -cov(W  ,W,)/var(Wl) ;  (2.3) 

m  m      m 


wm   s.d.(w  )        1/0 

-^  -   ^L=  d-r2)1/2,  (2.4) 


where 


W     s.d. (W  ) 
m         m 


^  ^     cov(W  ,W  ) 

r  =  corr(W  W  )  =    _w  "  m  .  (2.5) 

mm      o^   a~ 

m  m 


There  are  a  number  of  important  points  to  be  noted  about 
this  technique: 

(i)   It  allows  one  to  use  known  analytical  results  (such  as 
the  expected  value  of  the  limiting  waiting  time  in  an 


M/M/l   queue)  to  reduce  the  variance  of  a  simulation. 
Such  use  of  analytical  results  is  a  basic  principle  of 
simulation, 
(ii)   It  can  be  extended  to  non-linear  controls  or  to  multiple 
control  by  more  than  one  control  variable, 
(iii)   The  great  art  in  the  technique  lies  in  finding  a  control 
whose  mean  value  is  known  and  which  is  highly  correlated 
with  the  estimator  which  is  being  controlled.   Thus  (2.4) 
says  that   |r|   must  be  close  to  0.9  to  reduce  the  stand- 
ard deviation   a~    to  one  half  of   o~    and  this 

W  W 

m  m 

generally  is  equivalent  to  reducing  the  required  sample 
size  for  a  given  precision  by  a  quarter.   Controls  which 
are  that  highly  correlated  with  the  estimate  are  not 
easy  to  find, 
(iv)   It  is  in  general  too  much  to  ask  that  the  correlation  and 
variances  in  (2.3),  (2.4),  and  (2.5)  be  known  analyti- 
cally.  They,  therefore,  must  be  estimated  from  the 
simulation  data  and  this  will  reduce  the  variance  reduc- 
tion which  is  attained. 

Now  using  this  method  to  control  the  regenerative  estimate 
given  in  the  introduction, 

W(n)  =  Y(n)/a(n)  , 

where   n   refers  to  a  fixed  number  of  cycles  in  the  queue,  one  runs 

into  the  same  problem  of  synchronization  as  with  the  antithetic 

case,  namely   n   cycles  in  the   M/G/l   queue  may  take  a  very  different 


number  of  waiting  times  to  achieve  than  are  needed  for  n   cycles 
in  the  M/M/l   queue.   Moreover,  the  correlation  between   Y'   and 
Y   will  be  weak  and  made  even  weaker  by  the  use  of  the  ratio 
estimate.   This  problem  becomes  rapidly  apparent  in  simulation 
studies  of  the  technique. 

The  following  technique  which  we  have  called  internal 

(within  block)  control  has  been  developed  to  overcome  this.   It 
is,  in  fact,  a  special  case  of  the  very  broad  technique  called 
concomitant  variables  in  Gaver  and  Thompson  (1973) ,  p.  5  88  and 
will  be  illustrated  only  for  point  estimation  of  E (W)   and  E(a). 
Its  extension  to  other  quantities  discussed  by  Crane  and  Iglehart 

(19  74a)  is  immediate  in  principle,  although  the  control  quantities 
discussed  below  may  be  different. 

2 .1   Internal  Control:   Basic  Ideas 

The  idea  of  an  internal  control  variable  is  simple.  In  the 
estimate  W(n) ,  the  averages  Y(n)  and  a(n)  contain  n  random 
variables   Y    and   a,   which  are  each  functions  only  of  the   a, 

+-Vi 

interarrival  and  service  times  occurring  in  the   k —  cycle  (or  b.p.) 

and  are  independent  of  the  other  interarrival  times  and  service  times 

Thus,  it  is  natural  to  use  some  function  of  these   2a,   random 

variables  to  control  each   Y.   and   a,  . 

k        k 

The  naive  application  of  this  idea  is  that  if  we  can  reduce 
the  variance  of  both  the  numerator  and  denominator  Y(n)   and   a(n) 
we  will  reduce  the  variance  of  W(n) ,   but  we  will  show  that  the 
situation  is  more  complex  than  this.   We  will  denote  a  function  of 
the  random  variables  in  the   k —  cycle  by   C(k)   but  will  generally 


drop  the  index.   In  general  we  also  use   C  (k)   to  denote  a  control 

for  the  numerator  (top)  in  the  ratio  estimator  and   C_,  (k)   to 

o 

denote  a  control  for  the  denominator  (bottom) . 

Typically,   C(k),   or  simply   C,   could  be  the  difference 

between  the  service  time   vD    and  the  time  to  arrival  of  the  next 

3k 
customer  from  the  arrival  of  the   3vth   customer,  namely   yfi 

This  difference  has,  in  a   GI/G/1   queue,  a  known  mean 
y   -  A     and  large  positive  values  of  this  function   C(k)   corre- 
spond to  large  values  of  Y,   and   a,  ,   and  vice  versa.   We  return 
to  specific  control  variables  and  their  computation  in  the  next 
sub-section . 

Note  that  one  can  control  either  the  top  or  bottom  of  the 
ratio  estimator,  or  both;  to  fix  ideas  assume  we  control  the  top 

and  have,  in  general,  an  internally  controlled  estimate 

n 


FT  I      <Yk  +  3T(CT-E(cT))} 


WCT(n)  =  — iS-J: ,  (2.6) 

a 

where   $  ,   as  in  the  usual  control  estimation  technique,  is  fixed 
so  as  to  minimize   var (W   (n) ) .   In  practice  it  is  usually  esti- 
mated from  the  simulation  data. 

Now  the  quantity   a2/E2{a},   where   a2  =  E{Z2}   is  given  at 
(1.2),  is  just  the  leading  term  in  the  asymptotic  expansion  for  the 
variance  of  the  ratio  estimator  W(n) ,   and  carries  over  to  the 
more  complicated  situation  (2.6)  to  give  (asymptotically  as   n  -*■  °°) 

2 

n  var(WCT(n))  -  a2  =  {fj^f}  (var (Y' +6TC^-a ' ) 


/E(Y))' 
(E(a)  f 


varUy-a')  +  3^,0^},      (2.7) 
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where   Y'=Y/E{Y},   a'=a/E{a},   and   C.J,  =  CT/E{Y} .   For  a 
derivation  of  this  result  see  Cramer  (1946),  p.  354,  eq.  (27.7.3). 

It  follows  from  (2.7),  as  before,  that  to  minimize  the 
asymptotic  variance  of  W   (n)   we  must  take 

cov(Y,-a'  ,C,p    cov(Y,CT)   E  (y)  cov(a,CT) 
~3T  =    vaHcp     =   var(CT)  "  E~W  var(CT)   *     (2*8) 

But  most  importantly  we  notice  that   C*   must  be  highly  correlated 
with  the  difference   Y'  -a*.   To  achieve  this  is  much  more  difficult 
than  finding  a  quantity  which  is  highly  correlated  with  either  of 
Y   or   a,   simply  because   Y   and   a   are  highly  correlated  and 
increase  together.   In  particular  if   a=  1,   which  has  a  high 
probability  if   p,   the  traffic  intensity  is  small,  then   Y  =  0. 

Note,  too,  that  E  (W)  =  E(Y)/E(a),   the  quantity  we  are  trying 
to  estimate  in  the  simulation,  appears  in  (2.8)  for  the  top  control 
regression  coefficient.   Also  similar  equations  to  (2.7)  and  (2.8) 
pertain  to  the  case  where  the  bottom  of  the  ratio  is  controlled 
(in  this  case   a) ,   and  simultaneous  equations  can  be  derived  for 
3m   and   $_.   if  both  top  and  bottom  are  controlled.   The  additional 
complexity  in  estimating   8m   and   3D   does  not  seem  to  be  justi- 
fied  by  simulation  results  (discussed  later)  which  also  show  that 
if  only  one  control  is  used  there  seems  little  to  choose  in  terms 
of  reduction  achieved  between  putting  it  on  the  top  of  the  bottom. 
In  both  cases  the  control  must  be  highly  correlated  with  the 
difference  between   Y   and   a. 
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The  above  results  are  generally  applicable  for  any  regenerative 
process  in  which  ratio  estimates  are  used.   The  choice  of   C   is, 
of  course,  the  art  in  the  design  of  a  simulation  with  variance 
reduction  and  is  considered  for  the  GI/G/1   and  specifically  the 
M/M/l   queue  in  the  next  section.   In  all  cases  this  choice  is 
limited  by  one1 s  ability  to  compute  analytically   E(C) . 

2 .2   Internal  Control:   Design  Considerations 

We  discuss  here  the  design  of  internal  controls  for  the 
M/M/l   queue,  the  ideas  being  applicable  to  the   GI/G/1   case  with 
the  proviso  that  the  computations  might  be  considerably  more  diffi- 
cult.  Simple  computations  of   E(C)   are  given  here  and  more  diffi- 
cult ones  in  the  Appendix;  we  do  not  distinguish  between  bottom 
and  top  controls,  since  both  must  be  correlated  with  the  difference 
Y  -  a,   and  we  drop  consideration  of  cycle  number,  since  all  variables 
are  within  the  cycles  which  have  identical  structure. 

Again,  we  are  considering  estimation  of  E (W) ,   but  of  the 
many  possible  controls,  those  listed  below  would  probably  work  as 
well  with  other  functions  of  W,   e.g.  percentiles.   The  controls 
are  listed  roughly  in  order  of  complexity  of  computation  and  of 
supposed  correlation  with  Y-a.   This  can  usually  only  be  guessed 
at  and  generally  the  more  elaborate  controls  which  might  have 
greater  correlation  with  Y-a   are  more  difficult  computationally. 

Superscripts  on   C   are  labels  to  differentiate  the  controls. 
We  have  discussed  above  the  difference   X,  =  vn  -  u, /   whose  moments 
are  simple  to  compute.   Then  we  have 
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C(1)  =  X*  =  W  =  0    if   a  =  1  (2.9) 

=  xl   if   a  >  2.  (2.10) 

It  is  easy  to  compute   E(C    )   for  M/G/l   queues  and 
possible  for  the   GI/G/1   queue.   Thus,  for  the   M/M/l   case  we 
have,  using  the  Markov  property  of  the  exponential  distribution, 


P{v0<u1>  =  P{a=l>  = 


[1-F  (y)]dF  (y) 
0 


MJ 


(2.11) 

which  goes  from  1  to  1/2  as   p  =  A/y   goes  from  0  to  1;  furthermore, 


P{a>  2}  =  1-jL.jL  .  (2.12) 

Now  given  that   X..  =  vn  -  u.   is  greater  than  zero,  the 
excess   v  -  u,   is  distributed  as  an  exponential  random  variable 
with  parameter   u.   Therefore 

E(X+,  =E(C(1>)  MXjL^ia^)  .       (2.13) 

The  variance  of   C     can  also  be  computed. 

One  would  generally  like  to  obtain  more  correlation  of   C 
and  Y-a  when   Y   is  large,  and  one  feels  this  can  be  done  by 
bringing  in  the  second  waiting  time.   Thus  we  have  as  control 
candidates 
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ox   =0,  if   a  =  1 

CU)  (2.14) 

=  Xx  +  X2,  if   a  >  2; 

C(3)   =  (X^+X2)+  =  W2;  (2.15) 


0,  if   a  =  1 

CT4;  =  (2.16) 


Xl+  <Xi+x2)+' 
(2) 


if   a  >  2. 


The  control   C     is   WQ  =  0   if   a  =  l,   it  is   W Q  +  W.,   if 
a  =  2,   and  it  is   W~  +  W.  +  X?   if   a>  3.   It  is  an  attempt  to  capture 
the  effect  of  the  first  two  waiting  times  without  the  additional 
computational  complexities  involved  in  computing  the  expectations 
of  C(3)   and  C(4). 

Simulation  results  show  that  C     and  C     give  very 

(2) 
little  more  control  than   C     for  which,  in  the   M/M/l   case, 


we  have 


E(C(2))  =  0 +  E(X^+X2"/a>2) 


=  I£_{E(x;+x;|xJ>0} 


using  (2.12)  and  the  fact  that   X2   is  independent  of   X  .   We  use 

the  notation   E{x,A>   for  E{X1  },   where   X   is  a  random  variable, 

A   an  event,  and   1_   the  indicator  function  of  A. 

A 

Similar  computations  go  through  for  C  and  C  ;  from 
these  we  will  need  later  the  following  illustrative  result  (M/M/l 
case)  . 
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P{a=2}  =  P{X^>0,X^+v1<u2> 

=  P{X^>0/vQ-u1+v;L<u2} 


2 

l+pUiVU2-AJ       I+p     Ty+TT        TI+pI 


=  ^{P(u„>Y}  =  ^.  ,„^n2  =    P     ,       (2.18) 


where   Y,   the  sum  of  two  exponential (y)   random  variables,  has 
a  Gamma (y, 2)   distribution. 

From  this  we  also  get  that 

P(n>l\    -  1  -  1       P    -  P2(2+P)  ,2  19, 

P{a>3>  -  1-j+p   (i+p)  J  "  (l+p)'   '  (2-19) 

Now,  none  of  these  four  controls  is  specifically  designed 
to  be  correlated  with  the  difference   Y  -  a.   As  a  result  they  work 
well  for   y   and   A   such  that  Y   takes  on  values  much  larger 
than   a. 

To  fix  this  one  might  try 

c(5)  _  c(i)  _  ^  ±   =    lf2f3   or  4,  (2.20) 

since  the  mean   E(C    )   is  easily  calculated  if  E (C    )   is 
known  for   i  =  1,2,3,  or  4,   and  E(a)   is  known.   But  E(a)   is 
l/(l-p)   for  any  M/G/l   queue  (Cohen,  1969) ;  approximations  for 
the  GI/M/1   case  are  discussed  in  the  Appendix. 

There  is  an  additional  problem  of  dimensional  stability 
involved  in  using   C*   ,   as  with   C*   ,  C*   ,  C    ,  C*   ,   in  that 
E(C    )   depends  on  both   y   and   p,   while   E  (a)   depends  only  on 
p.   Thus  control  is  not  uniform  across  the  whole  range  of  param- 
eter values. 


15 


To  avoid  this  one  use 

c(6)  m   c(i)  _  a/^  (2.21) 

This,  with   i  =  2,   was  found,  in  simulation  studies,  to  be 
the  most  successful  control  variable  in  that  it  obtained  a  variance 
reduction  which  was  uniform  in   y   and   p   (or   y   and  X)       and  its 
mean  value  is  fairly  simple  to  compute.   These  simulation  results 
are  discussed  in  the  next  subsection. 

Note  that  multiple  control  variables  using  any  of  the  above 


controls  can  be  used.   In  particular,  one  need  not  take  the  differ- 
ed 
(1) 


(2) 
ence  of  say,   C     and   a/y,   but  may  use  a  multiple  control. 


However,  the  fact  that  two  regression  coefficients,  say   3T 

(2) 

and   3T    must  be  estimated  from  the  simulation  data  makes  the 

possible  gain  in  variance  reduction  of  dubious  value. 

Note  also  that  since  all  the  control  comes  from  within 
cycles,  there  is  no  reason  that  the  confidence  interval  estimation 
techniques  referred  to  in  the  introduction  would  not  go  through 
for  the  variance-reduced  estimates.   This  has  been  verified  in 
simulation  runs. 

2 . 3   Internal  Control:   Simulation  Results 

It  is  not  possible  to  verify  analytically  what  variance 
reduction  will  be  obtained  via  the  several  internal  controls 
listed  in  the  previous  section,  or  to  get  an  idea  of  the  magnitude 
of  the  effect.   Even  for  something  as  simple  as   C     it  is  diffi- 
cult to  compute  analytically  the  correlation  between   C     and 
Y-a   for  the   M/M/l   queue,  and  this  is  what  is  required  in  the 
equation  (2.5)  to  find  the  variance  reduction. 
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Thus,  we  resorted  to  simulations  to  verify  the  amount  of 
variance  reduction  obtained  and  the  relative  effectiveness  of  the 
various  controls.   In  the  final  simulations  all  runs  were  performed 
on  IBM  System  360/67  computer  using  the  LLRANDOM  package  (Learmonth 
and  Lewis  (1973) )  which  generates  random  numbers  according  to  the 
scheme  given  by  Lewis,  Goodman,  and  Miller  (1969)  and  exponentially 
distributed  random  numbers  using  the  Marsaglia  "rectangle-wedge- 
tail"  method.   Tests  of  the  random  number  generator  are  given  in 
Learmonth  and  Lewis  (1974) . 

Of  the  extensive  simulation  checks  performed,  we  give  here 
only  a  summary  of  the  conclusions  and  one  detailed  tabulation  and 
one  short  tabulation  in  the  case  of  the  most  suitable  control. 

(1)  The  controls   C    ,  C     and   C     do  much  better  gener- 
ally than   C    ,   with  little  improvement  over   C     obtained 
by  use  of   C     and   C    .   We  say  generally  because  results 
vary  with  X      and   y   and  their  ratio   p. 

(2)  Subtracting  the  number  of  customers  served  in  a  busy 
period  generally  improves  the  variance  reduction.   By  making 
it  dimensionally  stable  as  in  (2.21)  with   i  =  2   we  obtain  a 
"variance  reduction"  measured  in  terms  of  ratios  of  standard 
deviations,  of  approximately  70%,  uniformly  over   A   and   y. 

This  is  roughly  equivalent  to  halving  the  number  of  cycles  (b.p.'s) 
that  one  must  simulate;  (0.7) 2~  .5.   Much  better  reductions  can 
be  obtained  for  smaller   p   (i.e.   p  =  0.25)  by  specially  designed 
controls;  the  point  is  that   C  '      using   C     works  even  out 
at   p  =  0.99,   where  variance  reduction  is  extremely  important. 
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Table  1  shows  results  obtained  by  simulating  an   M/M/l   queue 
out  to  n  =  2000  cycles  and  replicating  the  simulation  250  times  to 
estimate  the  variance  of  the  estimates   W(n),   W   (n)  ,   W   (n)  , 
where  we  drop  the   n   for  convenience.   Here,  we  have  specifi- 
cally that 

i  n 

k    I     Yk 
W^(n)  = 


i  I      [ak+3B(C(6)-E{C(6)»] 
n  k=l 


where 


C(6)  =  c(2)  -a/u.  (2.22) 


The  estimated  precision  (standard  deviations)  of  the  estimates  of 
E (W)   are  given  in  brackets  under  the  estimates. 

The  results  in  Table  1  are  for   p  =  0.5   and  three  values  of 
y;   the  results  are  not  very  different  at  different  values  of   p.   The 
case   p  =  0.99  is  given  in  Table  2.   The  second,  third  and  fourth 
columns  in  the  Tables  give  correlations  between  the  control  and 
Y-a   etc.,  from  which  the  theoretical  variance  reduction  can  be 
computed.   They  are  very  close  to  the  values  given  in  the  next 
to  last  column,  from  which  we  deduce  that  estimating   3T   and   3R 
affects  the  variance  reduction  only  slightly.   There  is  negligible 
effect  of  different  values  of   y   for  fixed   p  =  0.5   and  fixed 
P  =  0.99. 

Note  that  for  the  results  for   p=0.99  given  in  Table  2,  the 
variance  reduction  is  73%  (about  the  same  as  for   p  =  0.5).   For  the 
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case  where  the  control  is  on  the  top,  i.e.  for   Y,   the  variance 
reduction  is  not  quite  as  good  for  control  of   a.   Note  too  that 
the  estimated  values  appear  in  some  cases  to  be  at  least  three  or 
four  standard  deviations  from  the  true  mean.   This  is  because  the 
estimates   W,   W   and  W   can  be  seen  from  the  100  replications 
to  be  non-normal.   In  other  words,  for  high   p(0.99),   the  simula- 
tion needs  to  be  taken  out  further  than  2000  cycles. 
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3.   Internal  Stratified  Sampling 

Another  technique  for  variance  reduction  which  can  be 
potentially  used  with  the  regenerative  method  is  stratified  sampling 
For  a  brief  description  see  Kleijnen  (1974),  p.  110.   In  essence 
this  uses  analytical  information  in  the  following  way. 

If  we  can  stratify  or  partition  a  random  variable   X  by 
its  values  (or  those  of  a  concomitant  variable)  into   K   strata 

labeled  k  =  l,2, ,K ,   we  can  write  the  mean  of   X   and  the  sample 

mean   X   as,  respectively, 

K 
y  =  E(X)  =  I      PE(X  |  X e  strata  k) ;  (3.1) 

k=l   K 


X 


k      *      Xi  =  £(   *   Xi+-.-+   I      X.)         (3.2) 
n  i=l   1    n\str  1  x       str  K  V 

n.  X.  n   X 

-     I    ifir+--+     I     if  ir  >  (3.3) 

str  1  n   1        strK    nK 
where   n,  =  number  of   X. 's   observed  in  strata,  and  p,   is  the 

k  1  K 

probability  of  being  in  strata. 

Now,  if  the   p,  ' s   are  known  and  we  substitute  them  in  (3.3) 
for  n ,/n,   we  get   X   ,   a  stratified  estimate.   It  will  be  biased, 

X.  S  L 

since  the  divisions  of  the  sums  in  the  populations  are  random;  if 

the  numbers  observed  in  each  population  are  controlled  and  taken 

to  be   nplfnp2,  etc.,  we  have  what  is  called  a  proportionally 

sampled  estimate   X    with 

ps 

K 

var(X)  =  var(X)  -  J   p,  (y,-y)2/n,  (3.4) 

PS  k=1   K   k 
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where   y,  =E{X  |  X  e  strata  k}.   Thus,  because  of  the  use  of  prior 
analytic  information  we  always  get  variance  reduction. 

The  variance  of   X    is  not  analytically  tractable,  but 
early  studies  reported  in  Kleijnen  indicate  it  is  close  to  (3.4) 
if   n   and  all  the   n  ' s   are  sufficiently  large.   Our  simulation 
studies  with  stratifying   a   have  confirmed  this;  because  of  the  ease 
of  computation  of   P{a=l} ,P{a=2} ,P{a>3} ,   as  in  section  2.3,  it  is 
natural  to  stratify   a.   Considerable  variance  reduction  is  obtained, 
especially  for  smaller  values  of   p,   the  traffic  intensity. 

Unfortunately,  when  the  quantity   a   is  stratified  in  the 
bottom  (denominator)  of  the  estimator  W(n)   very  little  overall 
variance  reduction  is  obtained  unless   p   is  small.   We  do  not 
understand  this  lack  of  variance  reduction  but  apparently  the 
correlation  between  the  stratified  version  of   a(n)   and  Y(n) 
is  reduced.   Analytical  studies  of  this  effect  are  very  difficult. 

It  is  also  possible  to  use  a  stratified  estimate  of  a 

—  (6) 

in  W(n)   and  to  control  the  top,   Y(n),   with  say   C     using 

(2) 
C    .   This  works  well  for  small   p,   but  increases  the  variance 

as   p   approaches   1.   Again  it  is  difficult  to  understand  this 

effect. 
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4 .   Conclusions  and  Summary 

We  have  been  able  to  obtain  a  worthwhile  variance  reduction 
using  internal  control  variables,  for  the  regenerative  estimate  of 
the  limiting  value  of  the  mean  waiting  time  in  an   M/M/l   queue. 
This  reduction  is  obtained  uniformly  over  all  parameter  values. 
It  is  fairly  certain  that  the  technique  will  work  well  with  any 
GI/G/1   queue  or  other  regenerative  stochastic  processes  or  systems. 
Internal  stratified  sampling  schemes,  however,  did  not  work  nearly 
as  well. 

The  techniques  can  be  extended  to  other  stable  stochastic 
systems,  such  as  the  Markov  chains  considered  in  Crane  and  Iglehart 
(19  74b) .   In  that  case  the  computation  of  the  mean  values  of  the 
controls  is  simpler  because  of  the  structure  of  the  Markov  chain. 

The  main  problem  in  applying  the  internal  control  variance 
reduction  technique  seems  to  lie  in  the  fact  that  the  estimator 
proposed  by  Crane  and  Iglehart  (19  74a)  involves  a  ratio  of  two 
random  variables,  and  these  are  difficult  to  work  with  in  general. 

An  alternative  which  will  be  considered  later  is  to  use  the 
existence  of  regeneration  points  more  specifically  to  obtain  vari- 
ance reduction  with  the  classical  estimator  W   given  at  (2.1). 

m 

One  advantage  which  the  regenerative  estimator  W(n)   has  over 

W   is  the  ease  of  obtaining  confidence  interval  estimates  or  esti- 
m  3 

mates  of  the  precision  of  W(n)   and  W(n) .   This  is  not  a  draw- 
back if  the  simulation  is  large  and  more  than  one  (say  ten  or 
twenty)  realizations  of  the  queue  are  obtained. 

To  fix  ideas  note  that  we  can  write  W   as 

ra 
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,   m 

W  =  -   >   W. 

m   m.i0   ] 

wNdn)  , 

=  ™'k=l  Yk  + YN(m)+ll  ' 

where  N (m)   is  the  number  of  completed  busy  periods  in  the  queue 
in   [0,m],   Y   as  before  is  the  sum  of  the  waiting  times  in  the 
k —  cycle,  and  Y '  .  .  .   the  sum  of  the  waiting  times  in  the  last, 
incomplete  cycle. 

Now  it  is  possible  to  apply  internal  controls  to  each   Y, 
in  the  sum.   Problems  arise  in  estimating  the  coefficients   $   in 
the  control  because  they  involve  a  random  sum  of  random  variables. 

But  it  is  much  easier  to  find  a  control   C   for  Y.   rather  than 

k 

the  difference   Y  -  a,  ,   and  also  it  is  still  possible  to  apply 
external  controls  as  well  as  internal  controls. 

These  ideas  will  be  followed  up  in  a  later  paper. 
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APPENDIX 

To  implement  either  the  internal  control  or  stratified 
sampling  techniques  certain  theoretical  parameters  associated 
with  the   GI/G/1   queue  are  required.   In  this  appendix  we  shall 
indicate  the  values  of  these  parameters  in  so  far  as  they  can  be 
calculated.   These  values  are  either  well-known  or  easily  calcu- 
lated.  For  a  reference  to  the  known  formulas  see  Cohen  (19  69)  . 

We  begin  with  E{a,},   the  expected  number  of  customers 

served  in  a  busy  period.   For  the  general   GI/G/1   queue  recall 

that  we  let   X  =  v   .  -  u   and  S  =  X,  +  . .  .  +  X  ,   for  n>  1, 
n   n-1   n        n    1        n  ' 

with   S_  =  0.   Then   a..  =inf{n>0:  S  50}.   The  general  expression 
for  E{a.)   is  given  by 

00 

E{a  }  =  exp{  I      n_1P[S  >0]}, 
n=l 

an  impossible  expression  to  evaluate  in  general.   Another  useful 

expression  for  E{a  }   is 

E{a  }  =  1/P{W=0},  (A.l) 

where   W  is  the  stationary  waiting  time.   In  the  special  case  of 
M/G/l,   however,  we  have 

ECc^}  =  (1+p)"1. 

Now  for  the  queue   GI/M/1   we  can  use  (A.l)  and  the  stationary 
distribution  of  the  embedded  Markov  chain  to  conclude  that 
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where   6   is  the  root  inside  the  unit  circle  of 

z  -  U{y (1-z) }  =  0 


-su. 


with   U(s)  =  E{e     },   Re  s>0,   and   u,   is  an  exponential  (y) 

r.v.   It  is  easy  to  check  that   6  =  p   for  M/M/l   queues.   Daley  (1975) 

has  recently  proposed  the  approximation  to   6   given  by 

S  =  a^l-p)2  +  2(l-b"1)p+  (2b"1-l)p2/ 

where   a  =  P{u=0},   E{u..}  =  l,   and  b  =  E{u2}.   This  approximation 
gives  good  results  in  a  number  of  examples  calculated  by  Daley  (1975) 
and  may  be  useful  for  the  purposes  of  this  paper. 

Next  we  turn  to  the  computation  of  P{a  =1}   and  P{a  =2}. 
For  the   GI/G/1   case  we  have 


and 


P{a.j=l}  =  PfS^O} 


P{a1=2}  =  P{S1>0/S2<0}, 


both  of  which  can  be  worked  out  with  a  little  effort.   For  the  M/M/l 
queue 


and 


P{a;L=l}    =    (1+p)    1 


P{ai=2}    =   p(l+p)    3. 


For  the  M/G/l  queue 


P{a;L=l}  =  V(X) 


Av, 


where   VQ)=E(e    °)  ,      and  for  the   GI/M/1   queue 
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P{a1=l}  =1-  u(y)  , 

where   U(s)   is  given  above.   For  the   M/E, /l   queue  and  the   E  /M/l 
queue  the  value   P{a  =2}   can  be  calculated  with  some  effort.   As 
these  expressions  are  cumbersome  they  shall  be  omitted. 

Next  we  turn  to  various  partial  expectations  which  are 
needed  for  internal  control 

E{S^+S2/a1>2}  =  E{S1/S1>0} + E{S2/S1>0} 
and 

E{a  ,a  >2}  =  E{a  }  -  P{a  =1} . 

In  the  special  case  of  the   M/M/l   queue, 

E{S1,S1>0}  =  p/{y(l+p) }, 

E{S+,Sl>0}  =  [2(I^)2+TI£l^]y-1/ 


and 


E{a1/a1>2}  =  2p/{ (1-p) (1+p) } 
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